{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Screening curve analysis\n",
    "\n",
    "Compute the long-term equilibrium power plant investment for a given load duration curve (1000-1000z for z $\\in$ [0,1]) and a given set of generator investment options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pypsa\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generator marginal (m) and capital (c) costs in EUR/MWh - numbers chosen for simple answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "generators = {\"coal\" : {\"m\" : 2, \"c\" : 15},\n",
    "              \"gas\" : {\"m\" : 12, \"c\": 10},\n",
    "              \"load-shedding\" : {\"m\" : 1012, \"c\" : 0}}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The screening curve intersections are at 0.01 and 0.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0,1,101)\n",
    "df = pd.DataFrame({key : pd.Series(item[\"c\"] + x*item[\"m\"],x) for key,item in generators.items()})\n",
    "df.plot(ylim=[0,50],title=\"Screening Curve\", figsize = (9,5))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = pypsa.Network()\n",
    "\n",
    "num_snapshots = 1001\n",
    "n.snapshots = np.linspace(0,1,num_snapshots)\n",
    "n.snapshot_weightings = n.snapshot_weightings/num_snapshots\n",
    "\n",
    "n.add(\"Bus\",name=\"bus\")\n",
    "\n",
    "n.add(\"Load\",name=\"load\",bus=\"bus\",\n",
    "      p_set=1000-1000*n.snapshots.values)\n",
    "\n",
    "for gen in generators:\n",
    "    n.add(\"Generator\",name=gen,bus=\"bus\",\n",
    "          p_nom_extendable=True,\n",
    "          marginal_cost=float(generators[gen][\"m\"]),\n",
    "          capital_cost=float(generators[gen][\"c\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.loads_t.p_set.plot.area(title=\"Load Duration Curve\", figsize = (9,5), ylabel=\"MW\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.lopf(solver_name=\"cbc\")\n",
    "n.objective"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The capacity is set by total electricity required.\n",
    "\n",
    "**NB:** No load shedding since all prices are below 10 000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.generators.p_nom_opt.round(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.buses_t.marginal_price.plot(title=\"Price Duration Curve\", figsize = (9,4))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prices correspond either to VOLL (1012) for first 0.01 or the marginal costs (12 for 0.49 and 2 for 0.5)\n",
    "\n",
    "**Except** for (infinitesimally small) points at the screening curve intersections, which correspond to changing the load duration near the intersection, so that capacity changes. This explains 7 = (12+10 - 15) (replacing coal with gas) and 22 = (12+10) (replacing load-shedding with gas). \n",
    "\n",
    "Note: What remains unclear is what is causing \\l = 0... it should be 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.buses_t.marginal_price.round(2).sum(axis=1).value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.generators_t.p.plot(ylim=[0,600],title=\"Generation Dispatch\", figsize = (9,5))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demonstrate zero-profit condition.\n",
    "\n",
    "1. The total cost is given by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(n.generators.p_nom_opt*n.generators.capital_cost \n",
    " + n.generators_t.p.multiply(n.snapshot_weightings.generators,axis=0).sum()*n.generators.marginal_cost)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. The total revenue by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(n.generators_t.p.multiply(n.snapshot_weightings.generators, axis=0)\n",
    " .multiply(n.buses_t.marginal_price[\"bus\"], axis=0).sum(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, take the capacities from the above long-term equilibrium, then disallow expansion.\n",
    "\n",
    "Show that the resulting market prices are identical.\n",
    "\n",
    "This holds in this example, but does NOT necessarily hold and breaks down in some circumstances (for example, when there is a lot of storage and inter-temporal shifting)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.generators.p_nom_extendable = False\n",
    "n.generators.p_nom = n.generators.p_nom_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.lopf();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.buses_t.marginal_price.plot(title=\"Price Duration Curve\", figsize = (9,5))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.buses_t.marginal_price.sum(axis=1).value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demonstrate zero-profit condition. Differences are due to singular times, see above, not a problem\n",
    "\n",
    "1. Total costs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(n.generators.p_nom * n.generators.capital_cost \n",
    " + n.generators_t.p.multiply(n.snapshot_weightings.generators, axis=0).sum() * n.generators.marginal_cost)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Total revenue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(n.generators_t.p.multiply(n.snapshot_weightings.generators, axis=0)\n",
    " .multiply(n.buses_t.marginal_price[\"bus\"],axis=0).sum())"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
